//    Copyright (C) Kherty.  All rights reserved.
using System;
using System.Collections.Generic;

namespace OpenLS.Spreadsheet.StandardAddIn.Statistics
{
    internal class Statistician
    {
        private const double DefaultEpsilon = 10e-9;
        private const double EPS = 2.22e-16;
        private const double LogGammaXMaxValue = 2.55e305;
        private const double Precision = 4.0*EPS;
        private const double XMININ = 2.23e-308;
        private static readonly double HALF_LOG_2_PI = 0.5*Math.Log(2.0*Math.PI);

        private static readonly double[] lanczos =
            {
                0.99999999999999709182,
                57.156235665862923517,
                -59.597960355475491248,
                14.136097974741747174,
                -0.49191381609762019978,
                .33994649984811888699e-4,
                .46523628927048575665e-4,
                -.98374475304879564677e-4,
                .15808870322491248884e-3,
                -.21026444172410488319e-3,
                .21743961811521264320e-3,
                -.16431810653676389022e-3,
                .84418223983852743293e-4,
                -.26190838401581408670e-4,
                .36899182659531622704e-5,
            };

        private static readonly int MaxIterations = 1000;

        internal static double MeanDifference(double[] sample1, double[] sample2)
        {
            return SumDifference(sample1, sample2)/sample1.Length;
        }

        internal static double SumDifference(double[] sample1, double[] sample2)
        {
            int n = sample1.Length;
            if (n != sample2.Length || n < 1)
            {
                throw new ArgumentException("array2.Length");
            }
            double result = 0;
            for (int i = 0; i < n; i++)
            {
                result += sample1[i] - sample2[i];
            }
            return result;
        }


        internal static double Mean(IEnumerable<double> numbers)
        {
            double result = 0;
            int count = 0;
            foreach (double d in numbers)
            {
                result += d;
                count++;
            }
            return result/count;
        }

        internal static double Covariance(double[] array1, double[] array2)
        {
            double average1 = Mean(array1);
            double average2 = Mean(array2);
            double result = 0;
            for (int i = 0; i < array1.Length; i++)
            {
                result += (array1[i] - average1)*(array2[i] - average2);
            }
            return result/array1.Length;
        }


        internal static double Variance(IEnumerable<double> numbers)
        {
            double average = Mean(numbers);
            double result = 0;
            int count = 0;
            foreach (double dbl in numbers)
            {
                double diff = dbl - average;
                result += diff*diff;
                count++;
            }
            return result/(count-1);
        }

        internal static double VarianceDifference(double[] sample1, double[] sample2,
                                                  double meanDifference)
        {
            double sum1 = 0d;
            double sum2 = 0d;
            int n = sample1.Length;
            if (n < 2 || n != sample2.Length)
            {
                throw new ArgumentException("sample2");
            }
            for (int i = 0; i < n; i++)
            {
                double diff = sample1[i] - sample2[i];
                sum1 += (diff - meanDifference)*(diff - meanDifference);
                sum2 += diff - meanDifference;
            }
            return (sum1 - (sum2*sum2/n))/(n - 1);
        }

        internal static double StandardDeviation(double[] array)
        {
            return Math.Sqrt(Variance(array));
        }

        public static double ErrorFunction(double x)
        {
            return Math.Abs(RegularizedGammaP(0.5, x*x, 1.0e-15, 10000));
        }


        public static double LogGamma(double x)
        {
            double result;

            if (double.IsNaN(x) || (x <= 0.0))
            {
                result = double.NaN;
            }
            else
            {
                double g = 607.0/128.0;

                double sum = 0.0;
                for (int i = lanczos.Length - 1; i > 0; --i)
                {
                    sum = sum + (lanczos[i]/(x + i));
                }
                sum = sum + lanczos[0];

                double tmp = x + g + .5;
                result = ((x + .5)*Math.Log(tmp)) - tmp +
                         HALF_LOG_2_PI + Math.Log(sum/x);
            }

            return result;
        }

        /**
		 * Returns the regularized gamma function P(a, x).
		 */

        public static double RegularizedGammaP(double a, double x)
        {
            return RegularizedGammaP(a, x, DefaultEpsilon, 1000000);
        }

        public static double RegularizedGammaP(double a,
                                               double x,
                                               double epsilon,
                                               int maxIterations)
        {
            double result;

            if (double.IsNaN(a) || double.IsNaN(x) || (a <= 0.0) || (x < 0.0))
            {
                result = double.NaN;
            }
            else if (x == 0.0)
            {
                result = 0.0;
            }
            else if (a >= 1.0 && x > a)
            {
                // use regularizedGammaQ because it should converge faster in this
                // case.
                result = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
            }
            else
            {
                // calculate series
                double n = 0.0; // current element index
                double an = 1.0/a; // n-th element in the series
                double sum = an; // partial sum
                while (Math.Abs(an) > epsilon && n < maxIterations)
                {
                    // compute next element in the series
                    n = n + 1.0;
                    an = an*(x/(a + n));

                    // update partial sum
                    sum = sum + an;
                }
                if (n >= maxIterations)
                    throw new ConvergenceException();
                else
                    result = Math.Exp(-x + (a*Math.Log(x)) - LogGamma(a))*sum;
            }

            return result;
        }

        public static double regularizedGammaQ(double a,
                                               double x,
                                               double epsilon,
                                               int maxIterations)
        {
            double result;

            if (double.IsNaN(a) || double.IsNaN(x) || (a <= 0.0) || (x < 0.0))
            {
                result = double.NaN;
            }
            else if (x == 0.0)
            {
                result = 1.0;
            }
            else if (x < a || a < 1.0)
            {
                // use regularizedGammaP because it should converge faster in this
                // case.
                result = 1.0 - RegularizedGammaP(a, x, epsilon, maxIterations);
            }
            else
            {
                // create continued fraction
                ContinuedFraction cf;
                cf = new ContinuedFraction(delegate(int n, double xx) { return ((2.0*n) + 1.0) - a + xx; }
                                           ,
                                           delegate(int n, double xx) { return n*(a - n); });

                result = 1.0/cf.Evaluate(x, epsilon, maxIterations);
                result = Math.Exp(-x + (a*Math.Log(x)) - LogGamma(a))*result;
            }

            return result;
        }

        public static double RegularizedBeta(double x, double a, double b)
        {
            return RegularizedBeta(x, a, b, DefaultEpsilon, int.MaxValue);
        }


        public static double RegularizedBeta(double x, double a, double b, double epsilon, int maxIterations)
        {
            double result;

            if (double.IsNaN(x) || double.IsNaN(a) || double.IsNaN(b) || (x < 0) ||
                (x > 1) || (a <= 0.0) || (b <= 0.0))
            {
                result = double.NaN;
            }
            else if (x > (a + 1.0)/(a + b + 2.0))
            {
                result = 1.0 - RegularizedBeta(1.0 - x, b, a, epsilon, maxIterations);
            }
            else
            {
                ContinuedFraction fraction;
                fraction = new ContinuedFraction(
                    delegate { return 1.0; },
                    delegate(int n, double xx)
                        {
                            double r;
                            double m;
                            if (n%2 == 0)
                            {
                                m = n/2.0;
                                r = (m*(b - m)*xx)/((a + (2*m) - 1)*(a + (2*m)));
                            }
                            else
                            {
                                m = (n - 1.0)/2.0;
                                r = -((a + m)*(a + b + m)*xx)/
                                    ((a + (2*m))*(a + (2*m) + 1.0));
                            }
                            return r;
                        }
                    );
                result = Math.Exp((a*Math.Log(x)) + (b*Math.Log(1.0 - x)) -
                                  Math.Log(a) - LogBeta(a, b, epsilon, maxIterations))*1.0/
                         fraction.Evaluate(x, epsilon, maxIterations);
            }
            return result;
        }

        private static double LogBeta(double a, double b, double epsilon, int maxIterations)
        {
            if (double.IsNaN(a) || double.IsNaN(b) || (a <= 0.0) || (b <= 0.0))
                return double.NaN;
            else
                return LogGamma(a) + LogGamma(b) - LogGamma(a + b);
        }

        public static double IncompleteBeta(double x, double p, double q)
        {
            if (x <= 0.0)
                return 0.0;
            else if (x >= 1.0)
                return 1.0;
            else if (p <= 0.0 || q <= 0.0 || (p + q) > LogGammaXMaxValue)
                return 0.0;
            else
            {
                double beta_gam = Math.Exp(-LogBeta(p, q) + p*Math.Log(x) + q*Math.Log(1.0 - x));
                if (x < (p + 1.0)/(p + q + 2.0))
                    return beta_gam*BetaFraction(x, p, q)/p;
                else
                    return 1.0 - (beta_gam*BetaFraction(1.0 - x, q, p)/q);
            }
        }

        private static double BetaFraction(double x, double p, double q)
        {
            int m;
            double sum_pq, p_plus, p_minus, c = 1.0, d, delta, h, frac;
            sum_pq = p + q;
            p_plus = p + 1.0;
            p_minus = p - 1.0;
            h = 1.0 - sum_pq*x/p_plus;
            if (Math.Abs(h) < XMININ)
                h = XMININ;
            h = 1.0/h;
            frac = h;
            m = 1;
            delta = 0.0;
            while (m <= MaxIterations && Math.Abs(delta - 1.0) > Precision)
            {
                int m2 = 2*m;
                // even index for d 
                d = m*(q - m)*x/((p_minus + m2)*(p + m2));
                h = 1.0 + d*h;
                if (Math.Abs(h) < XMININ)
                    h = XMININ;
                h = 1.0/h;
                c = 1.0 + d/c;
                if (Math.Abs(c) < XMININ)
                    c = XMININ;
                frac *= h*c;
                // odd index for d
                d = -(p + m)*(sum_pq + m)*x/((p + m2)*(p_plus + m2));
                h = 1.0 + d*h;
                if (Math.Abs(h) < XMININ)
                    h = XMININ;
                h = 1.0/h;
                c = 1.0 + d/c;
                if (Math.Abs(c) < XMININ)
                    c = XMININ;
                delta = h*c;
                frac *= delta;
                m++;
            }
            return frac;
        }


        internal static double LogBeta(double p, double q)
        {
            if (p <= 0.0 || q <= 0.0 || (p + q) > LogGammaXMaxValue)
                return 0;
            return LogGamma(p) + LogGamma(q) - LogGamma(p + q);
        }

        public static double BinomialCoefficientDouble(int n, int k)
        {
            return Math.Floor(Math.Exp(BinomialCoefficientLog(n, k)) + 0.5);
        }

        public static double BinomialCoefficientLog(int n, int k)
        {
            if (n < k)
                throw new ArgumentOutOfRangeException("n < k");
            if (n < 0)
                throw new ArgumentOutOfRangeException("n");
            if ((n == k) || (k == 0))
                return 0;
            if ((k == 1) || (k == n - 1))
                return Math.Log(n);
            double result = 0;

            // n!/k!
            for (int i = k + 1; i <= n; i++)
            {
                result += Math.Log(i);
            }

            // divide by (n-k)!
            for (int i = 2; i <= n - k; i++)
            {
                result -= Math.Log(i);
            }

            return result;
        }


        public static double FactorialDouble(int n)
        {
            if (n < 0)
            {
                throw new ArgumentOutOfRangeException("n");
            }
            return Math.Floor(Math.Exp(FactorialLog(n)) + 0.5);
        }

        public static double FactorialLog(int n)
        {
            if (n < 0)
                throw new ArgumentOutOfRangeException("n");
            double result = 0;
            for (int i = 2; i <= n; i++)
            {
                result += Math.Log(i);
            }
            return result;
        }
    }
}